function sac_text(varargin)

%fns = {'SP_an75.txt', 'SP_an64.txt', 'SP_an53.txt', ...
%    'SP_pyr75.txt', 'SP_pyr64.txt', 
fns=    {'SP_pyr53.txt', ...
  'SP_cw75.txt', 'SP_cw64.txt', 'SP_cw53.txt'};
fns = {'NetPYR-cw53-00000.txt', 'NetCW-cw53-00000.txt', ...
    'NetPYR-cw64-00000.txt', 'NetCW-cw64-00000.txt', ...
    'NetPYR-cw75-00000.txt', 'NetCW-cw75-00000.txt'};
base = '/users/pmanis/desktop/IgorStuff/';
base = '/users/pmanis/desktop/';
if(nargin == 0)
    iwhich = 0; % 0 = -53, 1 = -64, 2 = -75
    fn = fns{9};
    one_sac(base, fn);
else
    for i = 1:length(fns)
        one_sac(base, fns{i});
    end;
end;
return;

function one_sac(base, fn)

%X0=load([base fn], 'ASCII');
X0=[];
i=1;
hf = fopen([base fn], 'r');
w=0;
while(~feof(hf))
    l = fgetl(hf);
    if(length(l) > 0 & l(1) ~= '%' & l(1) ~= '"')
        X0(i,1:2) = sscanf(l, '%f %f');
        if(X0(i,2) ~= w)
            w = X0(i,2);
            i = i + 1;
            X0(i,1:2) = [NaN NaN];
        end;
        
            i = i + 1;
    end;
end;
fclose(hf);

b=find(isnan(X0(:,1)));
bs = 1;
dt=1; % 0.0125;
for j = 1:length(b)
    X{j} = X0(bs:b(j)-1,1)*dt;
    fprintf(1, '%d %9.2f\n', j, max(X{j}));
    bs = b(j) + 1; % skip the next NaN
end;

size(X)


twin = 100;
binw = 0.03333;
start1 = 00;
dur1 = 500;

[y, yh1, hx1, mr1] = sac(X, twin, binw, start1, dur1);
% perform gaussian fits on the histograms to identify the peaks and measure
% them. Using Molitor's routines (mrqfit).
% [...] = MRQFIT(F, P, X, Y, SIG, VP, LB, UB, IMAX, TOL)
% F = 'gaussian', P = [A0 A1 M1 S1 A2 M2 S2 ... AN MN SN]
% SIG is sigma Y (def = 1); VP is vary array per parameter [0 or 1];
% LB, UB are upper and lower bounds.
%
hx1c = hx1+0.5*(hx1(2)-hx1(1));
gpar1 = [1, 5, 0, 1]; % single gaussian centered on 0
scf = 1;
gulim1 = [0, 100, 10*scf, 20*scf];
gllim1 = [0, 0, 0.0, 0.01];
gvar1 =  [0 1 0 1];
nitermax = 400;

[fp11, chisq11, niter11, fitc11, err11, dep11] = mrqfit('gaussian', gpar1, hx1c*scf, yh1, [], gvar1, gllim1, gulim1, nitermax, []);



yg1 = gaussfunc(hx1c, fp11);
fwhmfac = 2*sqrt(2*log(2)); % note - log is ln.
fwhm1 = fwhmfac*fp11(4); % full width at half maximal height.


hf = findobj('tag', 'SAC2');
if(isempty(hf))
    hf = figure;
    set(hf, 'Tag', 'SAC2');
end;
% generate the results structure: SAC
SAC.hx1 = hx1; % save the histograms.
SAC.hy1 = yh1;
SAC.mr1 = mr1;

% gaussian fit results:
SAC.Gfit1 = fp11; % parameters
SAC.niter1 = niter11; % iterations
SAC.chisq1 = chisq11; % fit error
SAC.err1 = err11; % parameter estimate error
SAC.dep1 = dep11; % dependency between parameters
SAC.fwhm1 = fwhm1;
SAC.NPH1 = max(yg1);
%

figure(hf);
clf;


% text area
subplot('position', [0.1, 0.90, 0.8, 0.095]);
axis([0,1,0,1])
axis('off')
%ht(1)=text(0,0.80,sprintf('%-12s R[%d:%d]     %-8s  [%s]',DFILE.filename, DFILE.frec, DFILE.lrec, CONTROL(sf).protocol, date), 'Fontsize', 10);
%set(ht(1), 'interpreter', 'none'); % un-TeX the line - this is a filename and won't have tex chars, but might have an underscore.
%text(0,0.6,sprintf('Solution:%-12s  gain:%4.1f  LPF:%4.1f kHz', CONTROL(sf).solution, DFILE.igain, DFILE.low_pass(1)), 'FontSize', 8);
%text(0,0.4,sprintf('Ihold:%6.2f %s    RMP: %6.2f %s, Rin: %8.3f M\\Omega', ...
%    CONTROL(sf).iHold,CONTROL(sf).I_Unit, CONTROL(sf).Rmp, CONTROL(sf).V_Unit, CONTROL(sf).Rin), 'FontSize', 8);
text(0,0.200,sprintf('Window1: %.1f-%.1f  mean rate: %.2f s/s  ', start1, start1+dur1, mr1), 'FontSize', 8);
text(0, 0.000, sprintf('G1: A=%.2f (%.2f) S=%.3f (%.3f)  NPH = %.2f  FWHM1 = %.3f', ...
    fp11(2), err11(2), fp11(4), err11(4), SAC.NPH1, SAC.fwhm1), 'Fontsize', 8);


subplot('position', [0.1, 0.075, 0.8, 0.320]);
%bar(sqrt(hx1), yh1, 'histc');
semilogx(hx1c, yh1, 'ks', 'MarkerFaceColor', 'k', 'MarkerSize', 3.5);
hold on
semilogx(hx1c, yg1, 'r');
%semilogx(hx1c, fitc11, 'g');

u1 = get(gca, 'Ylim');
h1 = gca;
set(gca, 'Xlim', [0 100]);
xlabel('Delay (ms)');
ylabel('Normalized Peak Height');


y0 = u1(2);
x0 = 100;
axes(h1);
text(x0, y0, sprintf('Window 1 (%.1f - %.1f ms)', start1, start1+dur1), ...
    'horizontalalignment', 'right', 'verticalalignment', 'top', 'fontsize', 9);
box off


hout = fopen([base 'sac_' fn], 'w'); % output is input with sac_ prepended...
length(hx1c)
for i = 1:length(hx1c)
    fprintf(hout, '%f %f\n', hx1c(i), yh1(i));
end;
fclose(hout);

% now write the rasters...

    hout = fopen([base 'ras_' fn], 'w');
    t = ['ras_' fn];


done = 0;
row = 1;
while(~done)
    pflg = 0;
    if(row == 1)
        for i = 1:length(X)
            fprintf(hout, '"X%s_%d", "Y%s_%d", ', t, i, t, i);
        end;
        fprintf(hout, '\r\n');
    end;
    for i = 1:length(X)
        a=[X{i}]; % get the data itself
        if(length(a) >= row)
            fprintf(hout, '%f, %d, ', a(row), i);
            pflg = 1; % signal we printed at least one thing...
        else
            fprintf(hout, ', , '); % just the placeholders, no data
        end;
    end;
    fprintf(hout, '\r\n'); % end the row
    row = row + 1; % move to next row
    if(pflg== 0) done = 1; end;
end;


fclose(hout)




function [y] = gaussfunc(x, fp)
%
% calculate a gaussian based on x and FP
%
y = fp(1) + (fp(2)/(fp(4)*sqrt(2*pi)))*exp(-((x-fp(3)).^2)/(2*fp(4)^2));


function [y, yh, hx, mr] = sac(varargin)
%    [yan, sachist, sacx, mr] = sac(allsp, width, binw, delay, dur);
%
% Shuffled autocorrelation function
% Based on Louage et al., 2004
% X is an array of N responses x length(t) points in time. The times are in
% msec. Since the number of spikes in X may vary from trial to trial,
% X is assumed to be a cell array of spike times.
% The autocorrelation is calculated for every event in X that is not within
% is not within a single spike train.
% twin is the time window for the autocorrelation to be calculated (in
% msec)
% binw is the bin width for the histogram, in msec.
% delay is the delay to the start of measurement in X; in msec.
% dur is the duration of the measurement window in X; in msec.
% (these last two parameters allow you to pass a response and only analyze
% the portion that is relevant).
%
% y is the output autocorrelation function. It consists of a list of
% the correlation intervals, not corrected at all.
%
% yh is the histogram bin amplitudes, with bin width hx.
% This histogram is normalized.
% mr is the mean rate
%
% 10/5/04 Paul B. Manis, Ph.D.
% 10/5/04 some fixes for an spike trains.
% Note: calling the routine with letters as the first argument will generate the plots
% from Figure 2 of Louage et al.

global ALLCH DFILE


y = []; yh = []; hx = []; mr = 0;

if(nargin == 0) % assumes data is properly loaded already...
    return;
else
    X = varargin{1};
    twin = varargin{2};
    binw = varargin{3};
    delay = varargin{4};
    dur = varargin{5};
end;



if(nargin == 1) % call as a test - we make up our data and plug it in, and at the end
    % we make our own plot. the first arg is the letter corresponding to
    % the figures in Louage et al., 2004, Figure 2.
    %
    cmd = X; % capture command
    X={};
    baseper = 4/3;
    stimdur = 1000;
    ntestrep = 20;
    switch(cmd)
        case 'A'
            for i = 1:ntestrep
                X{i} = [baseper:baseper:stimdur];
            end;
            twin = 5;
            binw = 0.05;
            delay = 0;
            dur = stimdur;
            ddur = 10;
        case 'B'
            for i = 1:ntestrep
                X{i} = [baseper:baseper:stimdur];
                twin = 5;
                X{i} = X{i} + 0.08*randn(size(X{i}));
            end;
            binw = 0.05;
            delay = 0;
            dur = stimdur;
            ddur = 10;
        case 'C'
            fprintf(1, 'There is no computation for C, only A,B,D,E,F and G\n');
            return;

        case 'D'
            for i = 1:ntestrep
                X{i} = [baseper:baseper:stimdur];
                X{i} = X{i} + 0.170*randn(size(X{i}));
            end;
            twin = 5;
            binw = 0.05;
            delay = 0;
            dur = stimdur;
            ddur = 10;
        case 'E'
            for i = 1:ntestrep
                bl = [baseper:baseper:stimdur];
                bp = randperm(length(bl));
                X{i} = sort(bl(bp(1:floor(length(bp)/2)))); % elimnate half by random selection
                X{i} = X{i} + 0.170*randn(size(X{i}));
            end;
            twin = 5;
            binw = 0.05;
            delay = 0;
            dur = stimdur;
            ddur = 10;
        case 'F'
            twin = 5;
            binw = 0.15;
            delay = 0;
            dur = stimdur;
            for i = 1:ntestrep
                X{i} = stimdur*(sort(rand(120,1))); % corresponds to about 120 s/s
            end;
            ddur = 100;
        case 'G'
            twin = 5;
            binw = 0.15;
            delay = 0;
            dur = stimdur;
            sig = stimdur*(sort(rand(120,1)));
            for i = 1:ntestrep
                X{i} = sig; % corresponds to about 120 s/s
            end;
            ddur = 100;

        otherwise
            return;
    end;

end;

% now the SAC calculations.

lx = length(X);
% make sure input bin width and window are integer mulitples.
if(((twin/binw) - floor(twin/binw)) ~= 0)
    x=floor(twin/binw);
    twin = binw * x; % recalculate the window
end;
maxlx = 100;
maxn = 100000;
%fprintf(1, 'l = %d\n', lx);
if(lx > maxlx) % for testing purposes.
    lx = maxlx;
end;
yc=cell(lx,1); % Prellocation helps ALOT with speed
tic;
n = 1;
spcount = zeros(lx, 1);
for i = 1:lx
    xi = X{i}; % get spike train i
    kxi = find(xi >= delay & xi < delay+dur); % window the data to be analyzed
    xi = xi(kxi); % reduce data set
    spcount(i) = length(xi);
    fprintf(1, 'i=%d  ns: %d n', i, spcount(i)); % let user know we're working
    tic;
    yj = NaN*zeros(maxn,1); % preallocate space for cross correlation
    n = 0;
    for j = 1:lx
        if(j ~= i)
            xj = X{j}; % get spike train j
            kxj = find(xj >= delay & xj < delay+dur);
            xj = xj(kxj); % reduce data set
            for ti = 1:length(xi) % cross correlate against all spikes in xi
                iti = xi(ti); % time of ti'th spike in spike train i
                xjk = find(xj >= iti-binw/2 & xj < iti+twin+binw/2); % find spikes in j, in the window relative to current spike
                % binw/2 are egde effect corrections for the histogram stuff below.
                nx = length(xjk);
                if(~isempty(xjk) & (n+nx) < maxn)
                    n = n+1;
                    yj(n:n+nx-1) =  [xj(xjk)-iti]; % compute the difference in time between spike in j and relative to current spike in i
                    n = n + nx-1;
                end;
            end;
        end;
    end;
    yc{i} = yj; % somehow this is faster to store this way...
    fprintf(1, ' elapsed: %8.3fs  last n = %d\n', toc, n);
end;

y=[yc{:}]; % turn cell array into nx1 linear array
y=y(find(~isnan(y))); % clean it up - and shorten.

% now calculate the normalized histogram.
% normalization is N*(N-1)*r^2*deltat*D
% where N is the number of presentations (lx), r is the rate (sp/sec),
% deltat is the bin width for the histogram (sec), and D is the duration of
% the trace (?).
% normalization goes to 1 as the time gets larger...
%
rate = sum(spcount)/(lx*dur/1000); % get average rate
fprintf(1, 'Mean firing rate: %f\n', rate);
mr = rate;
nfac = lx*(lx-1)*rate*rate*(binw/1000)*(dur/1000); % correction factor
[yh] = histc(y, [0:binw:twin]); % generate histogram
hx = [0:binw:twin];
yh = yh/nfac; % to convert to rate, spikes/second
%m = mean(yh);
%fprintf(1, 'Mean value of yh: %8.3f, scale factor: %f\n', m, nfac);
if(nargin == 1) % if we are testing, generate a plot...
    hf = findobj('tag', 'sactest');
    if(isempty(hf) | (hf <= 0))
        hf = figure;
        set(hf, 'tag', 'sactest');
    end;
    figure(hf)
    clf;
    subplot(2,1,1);
    hold on;
    for j = 1:length(X)
        k = find(X{j} < ddur);
        plot(X{j}(k), j, 'ro', 'markersize', 2);
        set(gca, 'Xlim', [0 ddur]); % just show part of the train
    end;
    subplot(2,1,2)

    bar(hx, yh);
end;

return;


